***************************************************************************
* Set default directory with location of merged 2014 and 2019 data

cd "......\JMBE replication data\"
use 2019_2014_merged_MS_data.dta

svyset  [pw=final_weight], strat(pseudo_stfip)

* JMBE FIGURE 1
* Horizontal stacked bars
* Better labels
label define ad2 1 "Strongly agree"
label define ad2 2 "Agree", add
label define ad2 3 "Disagree", add
label define ad2 4 "Strongly disagree", add
label values emph_* ad2

catplot q1_globalwarm   [aw=final_weight], by(year, compact note("")) percent(year) recast(hbar) bar(1, bcolor(sienna)) blabel(bar, position(outside) size(4) format(%3.0f)) ytitle(" " "Percent teaching recent global warming")




* Generate plots without captions
* Figure 2a
catplot  emph_globaltemp year [aw=final_weight]  if q1_any_global==1 , /// 
aspectratio(.4) percent(year) scheme(s1mono) ///
bar(1, bcolor("53 33 0")) bar(2, bcolor(sienna)) /// 
bar(3, bcolor(sienna%80)) bar(4, bcolor(sienna%50)) /// 
title(" ", size(medium)) /// 
asyvars  stack var1opts(gap(0) sort((1)descending)) ///
var2opts(gap(*.2) label(labsize(medsmall)) ) outergap(*.2)  yla(0(10)100) /// 
legend(off) name(fig2a, replace) 

* Figure 2b
catplot  emph_sciconsensus year [aw=final_weight]   if q1_any_global==1 , /// 
aspectratio(.4) percent(year) scheme(s1mono) ///
bar(1, bcolor("53 33 0")) bar(2, bcolor(sienna)) /// 
bar(3, bcolor(sienna%80)) bar(4, bcolor(sienna%50)) /// 
title(" ", size(medium)) /// 
asyvars  stack var1opts(gap(0) sort((1)descending)) ///
var2opts(gap(*.2) label(labsize(medsmall)) ) outergap(*.2)  yla(0(10)100) /// 
legend(off) name(fig2b, replace) 
* Figure 2c
catplot  emph_naturalcauses year [aw=final_weight]   if q1_any_global==1 , /// 
aspectratio(.4) percent(year) scheme(s1mono) ///
bar(1, bcolor("53 33 0")) bar(2, bcolor(sienna)) /// 
bar(3, bcolor(sienna%80)) bar(4, bcolor(sienna%50)) /// 
title(" ", size(medium)) /// 
asyvars  stack var1opts(gap(0) sort((1)descending)) ///
var2opts(gap(*.2) label(labsize(medsmall)) ) outergap(*.2)  yla(0(10)100) /// 
legend(pos(6) col(2) size(medium)) name(fig2c, replace)

* Sig test of above:
svy, subpop(  if q1_any_global==1): tab emph_globaltemp year
svy, subpop(  if q1_any_global==1): tab emph_sciconsensus year
svy, subpop(  if q1_any_global==1): tab emph_naturalcauses year

** JMBE Figure 3
catplot typology_2016   if q1_any_global==1 [aw=final_weight], by(year, compact note("")) percent(year) recast(hbar) bar(1, bcolor(sienna)) blabel(bar, position(outside) size(4) format(%3.0f)) ytitle(" " "Percent")

*JMBE Table 1
balancetable wave ihave* if q1_any_global==1  [aw=final_w] using "balance_outcome_1.xlsx", replace vce(robust) pval ctitles("2014" "2019" "Difference" "p-value" "n") observationscolumn staraux wide(mean diff pval) format(%9.2f)



*JMBE Table 2 (Explaining change)

* A Standards --> Any hours
tabulate adopted_plus_2014 wave [aweight = final_weight]  , summarize(q1_any_global) nostandard noobs nofreq
svy: logit q1_any_global i.wave##i.adopted_plus_2014

* B Standard --> Mean hours (if any)
tabulate adopted_plus_2014 wave [aweight = final_weight]   if q1_any_global==1, summarize(q1_num_global) nostandard noobs nofreq
svy, subpop(  if q1_any_global==1): reg q1_num_global i.wave##i.adopted_plus_2014

* C Standards --> Emphasize warming
* reverse code
gen emp_global = (emph_globaltemp * -1)+5
tabulate adopted_plus_2014 wave [aweight = final_weight]   if q1_any_global==1, summarize(emp_global) noobs nofreq nostandard
svy, subpop(  if q1_any_global==1): reg emp_global i.wave##i.adopted_plus_2014

* D Standards --> Emphasize consensus
gen emp_consensus = (emph_sciconsensus * -1)+5
tabulate adopted_plus_2014 wave [aweight = final_weight]   if q1_any_global==1 , summarize(emp_consensus) noobs nofreq nostandard
svy, subpop(  if q1_any_global==1): reg emp_consensus i.wave##i.adopted_plus_2014

* Standards --> Equal time
tabulate adopted_plus_2014 wave [aweight = final_weight] if q1_any_global==1  , summarize(ihave_equaltime) noobs nofreq nostandard
svy, subpop( if q1_any_global==1): logit ihave_equaltime i.wave##i.adopted_plus_2014


* JMBE Table 3

* Demographics
gen white = 0
replace white = 1 if race_summary==1
balancetable wave female white seniority *_degree     [aw=final_w] using "balance_demo2.xlsx", replace vce(robust) pval ctitles("2014" "2019" "Difference" "p-value" "n") observationscolumn staraux wide(mean diff pval) format(%9.3f)

* Preparation
balancetable wave  prep_college prep_conted courseclimate_college sectionclimate_college courseclimate_conted sectionclimate_conted      [aw=final_w] using "balance_prep2.xlsx", replace vce(robust) pval ctitles("2014" "2019" "Difference" "p-value" "n") observationscolumn staraux wide(mean diff pval) format(%9.3f)

* Beliefs and attitudes
 gen believe_human_caused = .
 replace believe_human_caused = 1 if gw_causes==1
 replace believe_human_caused = 0 if gw_causes >1 & gw_causes < 99
 
 gen pid_dem = 0
replace pid_dem = 1 if partyid==1
replace pid_dem = . if partyid==.

gen pid_ind = 0
replace pid_ind = 1 if partyid==2
replace pid_ind = . if partyid==.

gen pid_rep = 0
replace pid_rep = 1 if partyid==2
replace pid_rep = . if partyid==.

gen bib_literal = 0
replace bib_literal = 1 if bible==1
replace bib_literal = . if bible==.

gen bib_inspired = 0
replace bib_inspired = 1 if bible==2
replace bib_inspired = . if bible==.

gen bib_ancientbook = 0
replace bib_ancientbook = 1 if bible==3
replace bib_ancientbook = . if bible==.

balancetable wave believe_human_caused consens80 pid_* govrole bib_*     [aw=final_w] using "balance_views2.xlsx", replace vce(robust) pval ctitles("2014" "2019" "Difference" "p-value" "n") observationscolumn staraux wide(mean diff pval) format(%9.3f)



* JMBE Table 4 - Mediation

* Table 4A - Preparation

mediate (q1_any_global adopted_plus_2014 , probit) (prep_college adopted_plus_2014, linear) (wave)  
estat proportion 

mediate (q1_num_global adopted_plus_2014 , linear) (prep_college adopted_plus_2014, linear) (wave)   if q1_any_global == 1
estat proportion 

mediate (emp_global adopted_plus_2014 , linear) (prep_college adopted_plus_2014, linear) (wave)   if q1_any_global == 1
estat proportion 

mediate (emp_consensus adopted_plus_2014 , linear) (prep_college adopted_plus_2014, linear) (wave)   if q1_any_global == 1
estat proportion 

mediate (ihave_equaltime adopted_plus_2014 , probit) (prep_college adopted_plus_2014, linear) (wave)  if q1_any_global == 1
estat proportion 

* Table 4B - Attitudes toward climate change
* Create warming attribution index
gen opinion = consens80 + gw_mostlyhuman

*For Linear Prob Model use mediate (q1_any_global, linear) ...
mediate (q1_any_global adopted_plus_2014 , probit) (opinion adopted_plus_2014 , linear) (wave)  
estat proportion 


mediate (q1_num_global adopted_plus_2014, linear) (opinion adopted_plus_2014, linear) (wave)   if q1_any_global==1
estat proportion 

mediate (emph_globaltemp adopted_plus_2014, linear) (opinion adopted_plus_2014, linear) (wave)    if q1_any_global==1
estat proportion 

mediate (emph_sci adopted_plus_2014, linear) (opinion adopted_plus_2014, linear) (wave)    if q1_any_global==1
estat proportion 

mediate (ihave_equaltime adopted_plus_2014 , probit) (opinion adopted_plus_2014 , linear) (wave)   if q1_any_global==1
estat proportion 